/*
 * DEM.h
 *
 *  Created on: 22/ago/2012
 *      Author: lordmzn
 */

#ifndef MOLEM_LANDSCAPEEVOLUTIONSTEP_H_
#define MOLEM_LANDSCAPEEVOLUTIONSTEP_H_

#include "DiscretizedSurface.h"
#include "Cell.h"
#include <vector>
#include <climits>
#include <unordered_map>

namespace MOLEM {
class UUID;
class DEMParam;

/**
 */
template<class T>
class LandscapeEvolutionStep {
public:
	/**
	 * Initialize LandscapeEvolutionStep with the given elevations of type T
	 */
	LandscapeEvolutionStep(std::vector<T>&, DEMParam&);
	~LandscapeEvolutionStep();
	/**
	 * check if DEM hasDepression=FALSE else calls depressionFilling,
	 * then calls the evaluateHydroNet to get the
	 * hydrological network of this DEM
	 */
	void evaluateHydroNet();

private:
	/* Pointer to const object: responsibility of demParam
	 * life is not taken by *this.
	 */
	DEMParam const* demParam;
	/**
	 * The DEM.
	 */
	Surface<SurfaceElement<T> > dem;
	/**
	 * maps cell into basins (the reverse is done by having
	 * a drainage basin field in Cell)
	 */
//	std::unordered_multimap<UUID, Cell<T> > drainageBasins;
	/**
	 * maps cell into links (the reverse is done by having
	 * a drainage basin field in Cell)
	 */
//	std::unordered_multimap<UUID, Cell<T> > drainageLinks;
	/**
	 * Mass of the terrain modelled. It's the sum of the
	 * elevation values contained in elevationData.
	 */
	double mass; // meters
	/**
	 * Value of the minimal positive difference of altitude
	 * for Planchon depression filling
	 */
	static const int epsilon = 1;
	/**
	 * Value of the huge number for transientSurface in Planchon
	 * depression filling
	 */
	static const T hugeNumber; // = INT_MAX;
	static const int maxDepth = 2000;

	/**
	 * Fills the depressions in the dem, evaluates its mass and
	 * its hydrological network.
	 */
	int evaluateHydroNet();
	/**
	 * Depression filling function according to Planchon. It calls
	 * DFPlanchonStage1 to initialize the transient surface, then
	 * calls DFPlanchonStage2 to fill up all the depressions.
	 */
	void depressionFilling();
	Surface<SurfaceElement<T> > DFPlanchonStage1();
	void DFPlanchonStage2(Surface<SurfaceElement<T> >&);
	int dryUpwardCell(SurfaceElement<T> c, Surface<SurfaceElement<T> >&, int);
	/**
	 * Evaluates the mass of this DEM.
	 */
	void evaluateMass();
	/**
	 * Find the direction of flow routing.
	 */
	void drainageDirectionAssignment();
	/**
	 * Evaluates the flow direction of the river that starts from
	 * the cell given as argument
	 */
	void GD8(SurfaceElement<T>&);
	/**
	 * labels each drainage basin with a different UUID
	 */
//	void drainageBasinLabelling();
	/**
	 * The function evaluates the drained area of the next cell. It
	 * is calculated as: number of cells whose flow is drained by
	 * the current cell * area of a single cell. Then it calls himself
	 * until a Sink or an already evaluated cell is found.
	 */
	void accumulateDrainage(SurfaceElement<T>);
	/**
	 * Labels each unbranched path on the flow routing network
	 * with a different UUID
	 */
//	void channelLinkRecovery();
};

template<class T>
LandscapeEvolutionStep<T>::LandscapeEvolutionStep(std::vector<T>& el, DEMParam& p) :
		demParam(&p), dem(el), /*drainageBasins(), drainageLinks(),*/mass(0.0) {
}

template<class T>
LandscapeEvolutionStep<T>::~LandscapeEvolutionStep() {
}

template<class T>
void LandscapeEvolutionStep<T>::evaluateHydroNet() {
	depressionFilling();
	drainageDirectionAssignment();
//	drainageBasinLabelling();
	SurfaceElement<T> current;
	// for each cell with isASink() == false
	if (current.isASource()) {
		accumulateDrainage(current);
	}
//	channelLinkRecovery();
}

template<class T>
void LandscapeEvolutionStep<T>::depressionFilling() {
	Surface<SurfaceElement<T> > transientSurface(DFPlanchonStage1());
	DFPlanchonStage2(transientSurface);
	dem = transientSurface;
}

template<class T>
Surface<SurfaceElement<T> > LandscapeEvolutionStep<T>::DFPlanchonStage1() {
	Surface<SurfaceElement<T> > transientSurface(dem.nr, dem.nc, hugeNumber);
	std::vector<SurfaceElement<T>& > trans_border;
	transientSurface.getBorder(trans_border);
	std::vector<SurfaceElement<T>& > dem_border;
	transientSurface.getBorder(dem_border);
	// copy the elevation of each cell in the border in digitalElevationModel
	// in the border of the transientSurface
	std::vector<SurfaceElement<T>& >::iterator it_trans = trans_border.begin();
	for (std::vector<SurfaceElement<T>& >::iterator it_dem=dem_border.begin();
			it_dem!=dem_border.end(); ++it_dem) {
//			it_trans.setElevation(it_dem.getElevation());
//			++it_trans;
	}
	return transientSurface;
}

template<class T>
int LandscapeEvolutionStep<T>::dryUpwardCell(SurfaceElement<T> current,
		Surface<SurfaceElement<T> >& transient, int depth) {
	depth++;
	if (depth > maxDepth) {
		return depth--;
	}
	// foreach neighbors of current
	// if elevation of neighbour cell in transient is = to hugeNumber
	// if cell.elevation in dem >= current.elevation in transient + epsilon
	// set transient.neighbour.elevation = dem.neighbour.elevation
	// depth = dryUpwardCell(neighbour, transient, depth)
	depth++;
	return depth;
}

template<class T>
void LandscapeEvolutionStep<T>::evaluateMass() {
	// sum all the elevations, transform them into mass and store it in mass
}

template<class T>
void LandscapeEvolutionStep<T>::drainageDirectionAssignment() {
	SurfaceElement<T>& top;
	do {
		// find top = the cell with max elevations among the ones without
		// a drainage direction already evaluated and not sea cell
		top.markAsSource();
		GD8(top);
	} while (1/**/);
}

template<class T>
void LandscapeEvolutionStep<T>::GD8(SurfaceElement<T>& t) {
	SurfaceElement<T>& reference = t;
	SurfaceElement<T>& current = t;
	SurfaceElement<T>& next;
	SurfaceElement<T>& next_primary;
	Direction primary_dir;
	Direction new_primary_dir;
	SecondDir secondary_dir;
	SecondDir new_secondary_dir;
	// for the first iteration, the LSD is chosen but the 2ndDir is evaluated
	// in order to evaluate the accrued error
	primary_dir = dem.findLocalSteepestDirection(current);
	secondary_dir = dem.findSecondaryDirection(current, primary_dir);
	current = dem.getNeighbour(reference, primary_dir);
	reference.setDrainsInto(current);
	// start the iteration toward a sink or an already evaluated cell
	while (!current.isASink() && !current.isDrainageDirectionEvaluated()) {
		new_primary_dir = dem.findLocalSteepestDirection(current);
		// check if next.markAsSink() is required: should never be required since a sink
		// it's when a depression is reached (but at this point they should be filled) or
		// when a cell with elevation == 0 is reached (but this happen only after the
		// application of the control or during the depression filling)
		new_secondary_dir = dem.findSecondaryDirection(current, primary_dir);
		// choose if primary or secondary has to be choosen
		next_primary = dem.getNeighbour(current, primary_dir);
		if (primary_dir == new_primary_dir) {
			if (secondary_dir == new_secondary_dir) {
				// if the gradient between the secondary and the reference is > than
				// the one between the primary and the reference
				if (dem.getGradient(reference, dem.getNeighbour(current, secondary_dir))
						> dem.getGradient(reference, next_primary)) {
					next = dem.getNeighbour(current, secondary_dir);
				}
			}
		} else {
			next = next_primary;
		}
		// check if the reference should be updated
		if (new_secondary_dir == SecondDir::UPDATE_REFERENCE) {
			reference = next;
		}
		// update reference to current so we progress toward the sink
		current = next;
		primary_dir = new_primary_dir;
		secondary_dir = new_secondary_dir; // if == UPDATE... it's not a problem
	}
}

//template<class T>
//void DEM<T>::drainageBasinLabelling() {
//
//}

template<class T>
void LandscapeEvolutionStep<T>::accumulateDrainage(SurfaceElement<T> current) {
	SurfaceElement<T> next;
	next = dem.getNext(current);
	if (!next.isASink() && !next.isDrainageAccumulationEvaluated()) {
		next.setDrainageAccumulated(
				current.getDrainageAccumulated/* + param.areaCell*/);
		accumulateDrainage(next);
	}
}

//template<class T>
//void DEM<T>::channelLinkRecovery() {
//
//}
}
#endif /* MOLEM_LANDSCAPEEVOLUTIONSTEP_H_ */
